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We study energy transport in a chain of quantum harmonic and anharmonic oscillators where 
the anharmonicity is induced by interaction between local vibrational states of the chain. Using 
adiabatic elimination and numerical simulations with matrix product states, we show how strong 
interactions significantly slow down the relaxation dynamics (with the emergence of a new time 
scale) and can alter the properties of the steady state. We also show that steady state properties 
are completely different depending on the order in which the limits of infinite time and infinite 
interaction are taken. 

PACS numbers: 05.60.Gg, 05.30.Jp, 03.75.Lm 


Miniaturization of devices is bound to reach the limits 
set by quantum mechanics in terms of both device output 
and heat/energy management. Therefore characterizing 
and controlling energy transport in nanostructures is crit¬ 
ical for the future engineering and design of quantum 
materials and devices. Fundamental studies in both the¬ 
ory and experiments are required. A promising approach 
for these studies is to analyze prototypical models which 
can simulate condensed matter systems and can be imple¬ 
mented experimentally in a clean, controlled and highly 
tunable manner [J-Q . As recent theoretical proposals [j| 
and experimental achievements have shown, cold 

ions are a great candidate for studying energy flow. In 
fact the ions form a regular artificial crystalline structure 
and each ion can vibrate around its equilibrium position 
thus affecting the vibration of the neighboring ions via 
Coulomb interaction. The use of site-resolved addressing 
allows control of the probability distribution of the differ¬ 
ent local vibrational modes. Moreover, in these set-ups 
it is possible to measure, both locally and globally, the 
excitation of the vibrational modes. 


Energy and particle transport in quantum bound¬ 
ary driven one-dimensional chains have been studied for 
spins, fermions and non-interacting bosons. Some no¬ 
table effects include negative differential resistance [9r 
111, emergence of long range order [l3 |. many-body reso¬ 


nances related to extended states in the cavity array 0, 
existence of more steady states with quantum-number 
conserving baths Il3 l , study of mesoscopic fl5l] or time- 
dependent baths [161, 117] and exact solutions in systems 
whose bulk is integrable 18-23J. 


In this work we show the influence of strong interac¬ 
tions on energy flow, relaxation dynamics and on the 
steady state for a bosonic chain (as it can be realized 
in cold ions experiments iEl). We first analyze the 
scenario in which the site-specific tailored dissipation 
produced by simultaneous application of red and blue 
sideband lasers acts only at one extremity of the chain. 
Henceforth we refer to it as ‘tailored bath’ to differenti¬ 



FIG. 1: (Color online) Depiction of the system under study. 
A chain of oscillators in which the ones at the extremities 
are local harmonic oscillators while those in the middle are 
anharmonic. The local oscillators are described by bosonic 
modes characterized by a local potential pi and local interac¬ 
tion strength Ui. Vibrations in each oscillator can tunnel to 
neighboring sites via the tunneling Ji. The harmonic oscil¬ 
lators are in contact with an environment through couplings 
A+ and A“. 


ate it from a thermal dissipative bath. We show that in 
this first scenario the steady state is independent of the 
strength of the interaction applied to other sites, in con¬ 
trast to the case in which the system is connected to a 
realistic thermal bath. The scenario with one single site 
connected to the tailored bath is also ideal to show how 
strong interactions significantly slow down the relaxation 
dynamics generating a new time scale which increases 
with increasing interaction strength. We then proceed 
to study the boundary driven case (i.e. the sites at the 
boundaries are driven by two tailored baths). In this 
case we show how interactions break ballistic transport 
and can be used to tune (but not to eliminate) the lo¬ 
cal bosonic occupation as well as the particle and energy 
currents in the chain. We stress that, in the strongly in¬ 
teracting regime, the emergent non-equilibrium physics is 
very different from that of hard core bosons (i.e. bosonic 
particles with infinite repulsive interaction). 

The system we analyze is a chain of L cold ions forming 
a linear lattice where vibrational quanta are transported 
along the chain. Since the vibrational quanta of the local 
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ion oscillator are bosons, we will refer to them as vibron 
[ 4 ]. The dynamics can be described by 
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where p is the density matrix of the system, H is the 
Hamiltonian and T> is the dissipative part of the evolution 
Q • The Hamiltonian H is given by |2~il Infil l 
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where Ji is the tunneling amplitude from site l to site 
l ± 1, Ui is the local interaction amplitude (due to in¬ 
duced anharmonicity from external lasers 0) and pi is 
the local potential (while pi can be locally controlled, we 
keep it constant throughout this work pi = p). In the 
remainder of the article we will take Ji = J constant 
throughout the chain. The operators ai and respec¬ 
tively destroy or create a vibron at site l. The tailored 
baths act locally on the harmonic oscillators injecting and 
extracting vibrons from two extremities: 


v(p) = J2 


1=1, L 


E ( '2aipaj - {a]di,p }) 
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Here A+ (A ; ) is the heating (cooling) rate of the dissi¬ 
pation. The system that we study is schematically repre¬ 
sented in Fig jT] which shows, for example, two harmonic 
oscillators coupled to the tailored environment at the 
boundaries and a bulk of anharmonic oscillators which 
are generated by local modification of the potential. The 
size of reduced Hilbert space that we consider scales as 
S = d 2L where d is the number of local vibrational modes 
considered. In our simulations d > 10 which for L = 4 
corresponds to S > 10 s . We will thus use analytical 
and numerical methods based on adiabatic elimination 
and the time-dependent Matrix Product States 


(t-MPS) algorithm [31 34|. 

A single bath: The steady state, p ss , for a chain in 
which only one site at an extremity is connected to 
the tailored environment is found to be a product state 

Pss = <S >1 {j2nPn,i\ n ^)( rl ^\) where \n,l) identifies the 
n-th vibrational state on site l and p s r ^ l the probability of 


it being occupied |35|, |36|. The value of is indepen¬ 


dent of the site l , the local interaction Ui and the local 
potential pi. It only depends on the amplitudes of the 
tailored dissipative coupling and it is given by 


f, (A+/AD" 

’ l hu (! - Af/Ar )" 1 


(4) 


This is clearly shown in Fig [2] for the illustrative case of 
a two site chain in which the first site (from the left) 



FIG. 2: (Color online) Diagonal elements of the density ma¬ 
trix of the second site at steady state 2 - The continuous 
black line represents the theoretical prediction; numerical re¬ 
sults are represented by the symbols: black circle o for U — 0, 
red x for U = 3 J and green + for U = 6 J. The dissipation on 
the left sites are fiA)*" = 3 J, ftAj~ = 6 J. The dashed lines show 
the distribution for a local thermal distribution for U = 3J 
(green circles) and U = 6J (blue squares). 


is coupled to the tailored dissipation while the second 
site can have different values of the interaction. The fig¬ 
ure depicts the steady state occupation of the vibrational 
modes [37| on the right site for different values of the in¬ 
teraction U 2 = 0, 3, 6 respectively by the black circles o, 
red x and green +. We note that the numerical results 
agree exactly with the theoretical expectation without 
any fitting parameter. The dashed lines with filled circles 
(U 2 = 3) and squares (U 2 = 6) represent instead the val¬ 
ues of p* s 2 ^ th e final state was a local thermal state with 
the same “temperature” as the first site. In this case the 
occupation of the local vibrational levels depends on the 
energy gap between the different levels and thus on the 
strength of the local interaction (anharmonicity). Hence, 
the dissipation studied here is not a thermal bath but a 
tailored dissipative process which tries to force a certain 
probability distribution of occupations in the energy lev¬ 
els of the density matrix independent of their energy. 

At this point it is possible to highlight a feature of this 
tailored dissipative system which will also be very impor¬ 
tant in the following. Suppose that the tailored dissipa¬ 
tion would be driving the system towards a product of 
local thermal states with the same temperature. In that 
case, when the repulsive interaction becomes significantly 
large, the anharmonic oscillators could be very well de¬ 
scribed by just the bottom two levels (also known as the 
hard core bosons limit) because the population of any 
other level would be significantly suppressed. However 
the actual tailored dissipation used imposes a particular 
local distribution independent of the vibrational energy 
levels. This means that very high energy levels would 
still be highly occupied and the hard core bosons picture 
cannot be used. In other words, the limits lim^oo and 
limj/^oo do not commute. The fact that high energetic 
levels are occupied also strongly affects the relaxation 
dynamics of the system as we analyze in the following. 

Two sites dynamics: As initial condition we consider 
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FIG. 3: (Color online) Probability of occupying level 2 of the 
oscillator on site 2, (P 2 , 2 ) versus time. The parameters used 
are: U = 30 J ,hAf = 3J, ftAj" = 9 J for the pink stars; 
U = 40J, hA/ = 3 J, ftAj" = 9 J for the red filled squares; U = 
50 J, hA/ = 3J, fiA]~ = 9J for the purple circles; U = 50J, 
hA/ = 4J, fiA]~ = 12J for the green X; U = 50 J, hA/ = 5J, 
ftA)~ = 15J for the light-blue +. 


a chain prepared in a tensor product of local thermal 
states with Ui = 0 V l and with constant average den¬ 
sity (hi) = (h m ) = 1 V l,m. Here ( O) = tr(pO) is 
the trace over the density matrix p of the operator O. 
Subsequently, at time t = 0 the left site is coupled to 
a tailored dissipative process that would impose a lo¬ 
cal average density (h\(t = oo)) = 0.5 lower than the 
initial density. At the same time, the interaction on 
site 2, C/ 2 , is also abruptly switched on. In the regime 
J -C A) 1- , A /,l = 1 ,L, we can adiabatically eliminate 
the dynamics of the site in contact with tailored environ¬ 
ment, and only consider the time evolution of the rest of 
the system (which, in this case, is the second ion) 27- 


_30|,[35j. The equation of motion for the diagonal elements 


of the reduced density matrix on the second site, p n , 2 , is 


dp n , 2 (t ) 
dt 


( /d~ -I- A+ 

l Xn ,n —1 ' Xn+1 ,7 


Pn,2 (/) 


(5) 


A Xn+l,n Pn+l, 2 (t) T Xn,n — 1 Pn—1 ,2(/) 


where X l m, n = 2J 2 AJ r m/ U$n 2 + h 2 (A/ — A t ) . From 

Eq. dHI) , which can easily be solved numerically, we 
learn of the presence of a new time scale r = 
Ui/ [J 2 (A, - - A+)] (when U 2 > h (Af - A+)). This 
result is clearly confirmed by Figj3] in which we plot the 
probability of having two vibrons at the second site given 
by (P 2 , 2 ) for different values of both the interaction Ui 
and coupling to the environment A f. All the curves col¬ 
lapse on top of each other once the time is rescaled with 
the new, emerging, time scale r. 

Boundary driven case'. Finally, we consider the case in 
which the chain is coupled to two tailored environments 
at its extremities. Each boundary is coupled differently, 
such that vibrons and energy currents are imposed. We 
show that, unlike the scenario with only one site coupled 
to the tailored environment, the local interaction strongly 
affects the profile of the vibron occupation, the particle 
current and the energy current in the system. We observe 
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FIG. 4: (Color online) (a) n 2 and (b) local fluctuations R 2 as 
a function of U 2 \ the blue triangles are results from the nu¬ 
merical t-MPS method while the red circles are analytic pre¬ 
dictions by adiabatic elimination. The dissipation coefficients 
used on the first and third sites are TiA/ = 5J, hA 7 = 21.66J 
and hA J = 5J and hA% = 55J. The insets represent the 
values of the density or fluctuations minus their asymptotic 
values for U —* oo, hj 0 and K 2 0 . The blue dashed lines repre¬ 
sent the numerical results while the red continuous lines are 
power-law fit with the exponent — 2 . 


first the case of a chain composed of only 3 sites, hence 
the interaction is only applied to the central site, U 2 ^ 0 
while U\ = U 3 = 0. In Fig@Ja) and (b) respectively, we 
show how the local density h 2 = (h 2 ) ss and local fluctua¬ 
tions k 2 = (n|) ss — (h 2 )i s of the steady state ( ss ) change 
with the interaction strength U 2 . These plots show a re¬ 
markable agreement between the exact simulation of the 
3 sites (blue triangles) and the adiabatic elimination ap¬ 
proach (red circles). For 3 sites, the equations stemming 
from the adiabatic elimination are given by 

l+Xn+l,n PnA*) 

1 = 1,3 / 

.(7 = ±1 / 

Z—1,3 

+ XrCn-l Pn-lAt) ( 6 ) 

Z—1,3 


dp n , 2 (t) 

dt 


The particular structure of these equations make it easy 
to compute the steady state [35 ] 
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Pn, 2 


<2 n 
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E l ,-h 
Z=l,3 Xk+l,k 
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where pQ S 2 se t s normalization. This expression 
agrees with our numerical simulations performed with 
t-MPS [H]. The insets of FigU computed using Eq.®, 
show that the values of h 2 and R 2 approach their re¬ 
spective asymptotic value for U 2 —> oo, and 
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FIG. 5: (Color online) Local density (hi) ss (a) and local 
energy {hi) sa (b) versus the site number for a chain of 4 
sites. The values of the interaction strength U2 = U3 are 
172 = 0 red stars, U 2 = 5 J black circles, U 2 = 10J blue dia¬ 
monds, U 2 = 25J purple squares and U 2 = 35J brown trian¬ 
gles. The dashed line represent the results from an hardcore 
bosons study. Other parameters are hAf = J , fiAj" = 4.33J, 
hAs = J and fiAj = 11J. 


as a power law with an exponent which, as predicted 
from Eq.©, is equal to —2. The asymptotic values 
are rhf = /3i/[(l - / 32)(1 + /3i — P 2 )] « 0.215 and 
it 2 ° = n§° (1 + h) / (1 — # 2 ) — (^ 2°) 2 ~ 0.233, where 

fi — ^ 1 1 -U / ^ _1 Ti 

M L(Ar-A +) 2 ^ (A--A+) 2 J / L1A7-A +) 2 (A 3 --A+H’ 

and @2 = (A) 1 " + Alj") / (A]j~ + A 3 ). Note that a hardcore 

0.167 and itf c = 0. 


bosons approach would predict n 


HC 
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Using MPS we further investigate a chain of 4 ion sites 
which would allow us to study how the local vibron den¬ 
sity ( fii) ss and local energy ( hi) ss are modified from the 
ballistic transport (typical of Ui = 0 ) to a new regime in¬ 
duced by the local interaction. The local energy is defined 
as hi = J/2 (a+cq_i + a^a; + i + h.c.)+(Ui/2)hi(hi-l) + 
fihi. In Fig IS] we show that the interaction breaks the 
ballistic transport, the local occupation of the levels de¬ 
creases and moreover it is not constant along the chain. 
The local energy, at large Ui is dominated by the inter¬ 
action and can be larger than the non-interacting case 
Ui = 0. The dashed lines in Fig[5] show the result of a 
hardcore bosons analysis in which the limit U —1 00 is 
taken before t —> 00 . 


We then turn to study how the interaction changes 
the current flow in the system. We analyze both particle 
current j n and energy current j h . Their expression is 
derived using the continuity equation 


md l = j?. hl -jZ +1 + s u ( 8 ) 


where ff _ 1; is incoming current, jf l+1 outgoing (both 
related to the Hamiltonian dynamics) and S) is a source 


FIG. 6: (Color online) Particle current j n (blue circles) and 
energy current j h (red diamonds) versus the site number for a 
chain of 4 sites. The red dotted and the continuous blue lines 
represent the hardcore bosons limit respectively for the energy 
and the particle currents. Other parameters are ftA) 1 " = J, 
hA 7 = 4.33J, fiAg" = J, ftAj = 11J and number of sites 
L = 4. 


term (due to dissipation). This results in 


Ji,k =i ^{a\ai - a\a k ) 

3l-i,t = + Uihi)bi-ih\ 

- bibj-^Ui-ihi-i + Urn)) 


(9) 


_ Yn Jl ~ 2 ’ 1 + 


Ui 

2 h 


ji-u 


( 10 ) 


In steady state the current is independent of the site and 
we can use lighter notations j n and j h . Fig ] 6 ] shows how 
both currents are lowered by the interaction towards an 
asymptotic limit different from hardcore bosons predic¬ 
tions. 

Experimental realization: The system studied here can 
be implemented using a chain of cold ions in a linear 
trap [24]. Though the results obtained here are indepen¬ 
dent of the specificities of the ion involved, the actual 
values of J and Ui depend on the mass and internal en¬ 
ergy levels of the ions involved. For Barium ions, 138 Ba + , 
it is possible to achieve an on-site interaction strength 
Ui ~ 5J with reasonable laser power |25|. Although go¬ 
ing beyond 8 J may not be feasible u sing a standing wave 
configuration as proposed by Porras [24], the effects of in¬ 
teraction on transport is quite pronounced even at lower 
interaction strength (see Figs. [5] and [ 6 ]) • The tailored 
bath has been implemented in an ion trap setup [ 0 ] while 
hopping of transverse vibron has been implemented by 
Haze et al. [1]. Addressing individual sites for a chain 
of 4 ions with inter-ionic distances as large as tens of 
microns does not seem to pose any experimental chal¬ 
lenge, however for a large number of individually ad¬ 
dressed ions, laser beams with finite focussing may af¬ 
fect the neighbouring ions. Therefore we believe, with 
present day state-of-the-art techniques it is possible to 
realize the system proposed here as a proof-of-principle 
for a chain below ten ions. 
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Conclusions: We have shown that anharmonicity in¬ 
duced interactions break the ballistic transport of parti¬ 
cles and energy in the non-interacting regime and can be 
used to regulate the time scales and the steady state of 
energy transport in linear chains. We have developed an 
analytical understanding of these phenomena thanks to 
adiabatic elimination and shown that the properties of 
the system in the strong interaction limit are very differ¬ 
ent from a hardcore bosons approach. This work is an 
important step towards the understanding of energy and 
heat flow in interacting bosonic systems. 

We are thankful to U. Bissbort, S. Denisov, S. Flach, 
M. Fleischhaur, Z. Gang, N. Li, D. Rossini and C. Teo for 
insightful discussions. This project is funded by SUTD 
Start-up grant EPD2012-045. 
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SUPPLEMENTARY MATERIAL 


which is identical to Eq.(2S). 


Steady state for a chain with one boundary site 
under dissipation 


Adiabatic elimination 


In this section we aim to show that Eq. (4) of the main 
article is indeed a steady state of the system. First of all 
we note that V(p ss ) = 0. In fact 

V(pss) = E A 5" i “ 1 l 

“? Ar UXfMf n|n>(n| 

+ S A * i <A -‘i+/| r (» + i)l" +1>(" + il 

(1S) 

where the first line in Eq. (I1SI) cancels the fourth and sec¬ 
ond is equal and opposite to the third. It is also possible 
to show that [H,p ss \ = 0. The interacting and the po¬ 
tential part of the Hamiltonian are readily verified. As 
for the kinetic part of the Hamiltonian, it is sufficient to 
show that [ajd p , p ss ] = 0. Indeed you have 
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Here we show one way to derive Eq.(5) of the main 
paper. In the regime in which A 1 J it is possible to 
write classical diffusion equations for the evolution of the 
diagonal elements of the single site reduced density ma¬ 
trix p U: 2 - We will follow the method described in detail 
in the supplementary material of |3 qJ . Using the same 
notations as [30] we divide the Liouvillian evolution in 
two parts, 


C{p) = £ 0 (p) + V(p) 
i 

h . 
i 

h L 


4>(p) = - 

V(p) = - 


Hd,P 

Hk,P 


V(p) 


(4S) 

(5S) 

(6S) 


where H D = U 2 / 2 E„ 2 A 2 (n 2 - 1) + pE ni ,n 2 (™2 + ™i) 
and Hk = — J (b\b 2 + h.c.^ . Let us take 


P Q ° = (E^,2|n 2 )(n 2 | I 0pf 


= (Y,Pn,M(n 2 \j 0 (^E—l ni >< ni lj ( 7S ) 

where a = A^/A] - and /? = (l — A+ /A]”) 1 and p H2 ,i 
has the same meaning as p n p in the main paper, ff 10 
lives in the so-called Ho, decoherence free, space and it is 
connected to dissipating spaces Q a via the kinetic term 
of the Hamiltonian which is a small perturbation [39i ]. 
The time evolution of p n ° is given by 


dt 


p n ° = £ (p n °) (8S) 

£ = £o-^V(£ 0 ) _1 V (9S) 


where, in our case, £o (p Q °) = 0. V(p n °) gives 4 different 
terms and we choose one of them, C = ( iJ/h)b\b 2 p n °, 
to exemplify the derivation of Eq.(5) of the main paper. 
This gives 

C 0 (C) = (iU 2 (n 2 - 1) - (AZ - A+)) C (10S) 

This results in the diagonal terms of the density matrix 
to follow the differential equation 


= Af/Al 

0 \n) p (n + l\ 


_m,n=0 

(A+/A~) m+n 

(1-A+/Ar)- 2 J 


+ l)|m + 1 )z(to 



(3S) 

= (n 2 |tri^£p° 0 ^ |n 2 ) 

(US) 
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where tri traces out the first site. After some computa¬ 
tion this turns into 


d 

Jt Pn2 


,2 — 


2J 2 (n 2 + l)ni(A 1 — A+) 
t-l U 2 n 2 + h 2 ( A r -A+)2 


Since 


,n 1 — 1 


a 1 a 

x I Pn 2 ,2—5-Pn 2 + 1,2 


2 J 2 (m + l)n 2 (A 1 — A+) 
U 2 (ri 2 - i ) 2 + h 2 (Pq - A+) 2 

a 711 a ni+1 ' 

Pn 2 ,2~3 -Pn 2 —1,2 — 


P 



At 

Ar-A+ 




Ar 

Ar-A+ 


we obtain 

d 2J 2 (ri2 + l)A^“ 

df P ” 2 = “ f7 2 n 2 + ft 2 (A^ -A+) 2pn2 

2 J 2 n 2 A(” 

~ t/ 2 (n2-l) 2 + ft 2 (A^-A+) 2Pn2 
2J 2 (n 2 + 1)A+ 

+ f7 2 n 2 + ft 2 (Ar - A+)2 Pn2_1 
2 J 2 n 2 A(~ 

+ U 2 (ri 2 - l) 2 + ft 2 (A^ - A+) 2 Pn2+1 


(12S) 


(13S) 

(14S) 


(15S) 


which is Eq.(5) of the main paper. In an analogous man¬ 
ner Eq.(6) can be derived. 


Solving the first equation allows to compute 


Pi,2 = 


E Z,+ 

i=l,3 *1,0 ss 

^ — P 0,2 

2 ^, 1 — 1,3 Xi,o 


which substituted in Eq.(17S), gives 


ss _E/Z—1,3 *2,1 ss 

P 2,2 — — — Pi,2 

2^il= 1,3 *2,1 

(Ez= 1,3 *2,l) (Ez=1,3 *1,0 ) ss 

”(£z=i, 3 *2’,l) (Ez=i, 3 xi’,o ) A5,2 


(20S) 


(21S) 


Steady state for 3 sites 


Solving in sequence the following equations results in 


Solving Eq.(6) of the main paper for the steady state 
(which satisfies dp n , 2 ss /dt = 0) we get a set of coupled 
equations 

- E xi’>o,2 ss + e xi’,o pi :2 = 0 (i6S) 

1 = 1,3 1 = 1,3 


- E 

V=l,3 

+ *1,0 Po^2 = 0 

;=i,3 


0 + X2,i ] Pl?2 + E *2,1 P2% 


1 = 1,3 


(17S) 


E z,+ 

1=1,3 Xn,n—1 nSS 
Pn, 2 —— j~Z Pn- 1,2 

Ei=l,3 Xn,n— 1 


(Ez=1,3 Xn%- 1 ) ■ • ■ (Ez= 1,3 *l,o) 
(Ez=1,3 Xn,n- 1) • ■ • (Ez=1,3 *1,o) 


Po ,2 (22S) 


- E *E- 


1 + Xn+l ,n I Pn,2 


, 1=1,3 


(18S) 


E Xn+l,n Pn+ 1,2 + E Xn^n-l Pn- 1,2 = 0 
1= 1,3 1= 1,3 


which is Eq.(7) of the main paper. 


(19S) 



















